rm(list=ls())
library(pheatmap)
library(ggplot2)

'<-'(
  PCA,
  function(A, tag = NULL, center = T, normalize = F, d = 2, draw = T)
  {
    m <- dim(A)[1]
    n <- dim(A)[2]
    if(center)
    {
      A <- apply(A,2,function(x){return(x-mean(x))}) 
    }
    if(normalize)
    {
      A <- apply(A,2,function(x){return(scale(x))})
    }
    B <- 1/(m-1) * t(A) %*% A
    tmp <- eigen(B)
    values <- as.vector(tmp$values)
    disVaule <- values / sum(values)
    vectors <- as.matrix(tmp$vectors)
    vectors <- vectors[,1:d]
    C <- A %*% vectors
    if(d==2 && draw)
    {
      if(is.null(tag))
      {
        plot(C[,1],C[,2],xlab=paste('PC1(',disVaule[1],'%)'),ylab=paste('PC2(',disVaule[2]*100,'%)'))
      }
      else
      {
        types <- levels(factor(tag))
        classes <- lapply(types, function(x){return(which(tag==x))})
        plot(C[classes[[1]],1],
             C[classes[[1]],2],
             col=2,
             xlab=paste('PC1(',round(disVaule[1]*100,2),'%)'),
             ylab=paste('PC2(',round(disVaule[2]*100,2),'%)'),
             xlim=c(min(C[,1]),max(C[,1])),
             ylim=c(min(C[,2]),max(C[,2])),
             pch=16)
        classes <- classes[-1]
        col_cnt <- 2
        
        for(each in classes)
        {
          col_cnt <- col_cnt + 1
          lines(C[each,1],C[each,2],col=col_cnt,pch=16,type='p')
        }
        for(x in 1:dim(vectors)[1])
        {
          tname <- names(A[1,])[x]
          print(paste(c(getLetter(x),tname)))
          lines(vectors[x,1],vectors[x,2],col=1,pch=getLetter(x),type='p')
          lines(c(0,vectors[x,1]),c(0,vectors[x,2]),col=1,type='l')
        }
        legend('topleft',legend = types, col=2:col_cnt, pch=16)
      }
    } else {
      return(C)
    }
  }
)

getLetter <- function(id){
  wh = ceiling(id / 26 ) - 1
  etl = 'z'
  if (id %% 26 != 0){
    etl = letters[id %% 26]
  }
  if(wh == 0){
    return(etl)
  } else {
    return(paste(c(rep('z',wh),etl), collapse = ''))
  }
}

fromLetter <- function(id){
  id <- tolower(id)
  id <- strsplit(id, '')[[1]]
  if(length(id) <= 1){
    return(which(letters==id))
  } else {
    ix <- 0
    for(i in 1:length(id)){
      ix <- ix + which(letters==id[i]) * 26^(length(id)-i)
    }
    return(ix)
  }
}


'<-'(
  PCA,
  function(A, tag = NULL, center = T, normalize = F, d = 2, draw = T)
  {
    m <- dim(A)[1]
    n <- dim(A)[2]
    if(center)
    {
      A <- apply(A,2,function(x){return(x-mean(x))}) 
    }
    if(normalize)
    {
      A <- apply(A,2,function(x){return(scale(x))})
    }
    B <- 1/(m-1) * t(A) %*% A
    tmp <- eigen(B)
    values <- as.vector(tmp$values)
    disVaule <- values / sum(values)
    vectors <- as.matrix(tmp$vectors)
    vectors <- vectors[,1:d]
    C <- A %*% vectors
    if(d==2 && draw)
    {
      if(is.null(tag))
      {
        plot(C[,1],C[,2],xlab=paste('PC1(',disVaule[1],'%)'),ylab=paste('PC2(',disVaule[2]*100,'%)'))
      }
      else
      {
        types <- levels(factor(tag))
        classes <- lapply(types, function(x){return(which(tag==x))})
        plot(C[classes[[1]],1],
             C[classes[[1]],2],
             col=2,
             xlab=paste('PC1(',round(disVaule[1]*100,2),'%)'),
             ylab=paste('PC2(',round(disVaule[2]*100,2),'%)'),
             xlim=c(min(C[,1]),max(C[,1])),
             ylim=c(min(C[,2]),max(C[,2])),
             pch=16)
        classes <- classes[-1]
        col_cnt <- 2
        
        for(each in classes)
        {
          col_cnt <- col_cnt + 1
          lines(C[each,1],C[each,2],col=col_cnt,pch=16,type='p')
        }
        for(x in 1:dim(vectors)[1])
        {
          tname <- names(A[1,])[x]
          print(paste(c(getLetter(x),tname)))
          lines(vectors[x,1],vectors[x,2],col=1,pch=getLetter(x),type='p')
          lines(c(0,vectors[x,1]),c(0,vectors[x,2]),col=1,type='l')
        }
        legend('topleft',legend = types, col=2:col_cnt, pch=16)
      }
    } else {
      return(C)
    }
  }
)
'<-'(
  my_lm,
  function(x,y,d=1,draw=F,show=T,col=2)
  {
    ox <- x
    x <- sapply(0:d,function(v){return(x^v)})
    e <- try(w <- solve((t(x)%*%x))%*%(t(x)%*%y))
    if(typeof(e)=='character')
    {
      warning('not good dim')
      #return(my_lm(x,y,d-1,draw,show,col))
      return(F)
    }
    if(draw)
    {
      if(draw!=2)
      {
        plot(ox,y,col=3,xlab='x',ylab='y',main=paste('lm: dim=',d))  
      }
      tx <- seq(min(ox)-10,max(ox)*1.2,0.1)
      x <- sapply(0:d,function(v){return(tx^v)})
      lines(tx, as.vector(x%*%w), col=col, type='l')
      if(show){
        legend('topright',legend=c('y=Σki*x^i',w))
      }
    }
    return(w)
  }
)

'<-'(
  lm_t,
  function(x,y,w)
  {
    x <- cbind(1,x)
    yp <- x%*%w
    n <- dim(x)[1]
    p <- length(w)
    c <- solve(t(x)%*%x)
    tj <- numeric(p)
    s <- n-p-1
    t_sd <- sqrt((1/s)*sum((y-yp)^2))
    for(i in 1:p)
    {
      tj[i] <- -1 * w[i] / sqrt(c[i,i]*t_sd)
    }
    tj <- abs(tj)
    return(1-pt(tj,s))
  }
)

dat <- read.csv('./fenshuidiaoyan/wenjuan.csv')
ndat <- matrix(as.numeric(unlist(dat[,-c(1,81,82)])),134,79)

timu <- names(dat[1,])
pie_which <- c(2, 4, 5, 7, 8, 9, 10, 15, 17, 21, 22, 23, 25, 26, 27, 28)
my_col <- function(x){
  if(x > 6){
    return(1:x)
  }
  cols <- c(8,43,36,74,123,26)
  return(cols[sample(1:6, 2, replace = F)])
}


sapply(pie_which, function(each){
  pdat <- dat[,each]
  xuan <- levels(factor(pdat))
  weight <- sapply(xuan, function(x){length(which(pdat == x))})
  minzi <- (weight / sum(weight)) * 100
  minzi <- round(minzi, 2)
  minzi <- sapply(minzi, function(e){paste(as.character(e), '%',seq='')})
  new_col <- my_col(length(xuan))
  #pie(weight, labels = minzi, col=new_col, main=timu[each])
  if(length(xuan) == 2 && all(xuan == c('0','1'))){
    xuan <- c('否','是')
  }
  #legend('bottomright', legend=xuan, fill=new_col)
})

retu_which <- c('B','C','D','E','F','G','H','I','J','K','L','M','W','X','Y','AI','AO','AP','AR','AS','AT','AU','AV','BH')
retu <- sapply(retu_which, fromLetter)
cor_max <- sapply(retu, function(each){
  return(sapply(retu, function(another){
    return(cor(dat[,each], dat[,another]))
  }))
})
pheatmap(cor_max, display_numbers = matrix(ifelse(abs(cor_max < 0.05), "*", ""), nrow(cor_max)), cluster_row = FALSE, cluster_cols = F)
retu_name <- t(rbind(retu_which,timu[retu]))

infos_id <- c('B','D','E','F','G')
yins_id <- retu_which[-(1:6)]
infos <- sapply(infos_id, fromLetter)
yins <- sapply(yins_id, fromLetter)
pics <- list()
cnt <- 1
for(info in infos){
  aldat <- c()
  ceng <- levels(factor(dat[,info]))
  ceng_num <- sapply(ceng, function(t){length(which(dat[,info]==t))})
  for(yin in yins){
    if(length(levels(factor(dat[,yin]))) != 2){
      next
    }
    pdat <- dat[which(as.character(dat[,yin])=='1'),info]
    ceng <- levels(factor(pdat))
    weight <- sapply(ceng, function(t){length(which(pdat==t))/ceng_num[t]})
    aldat <- rbind(cbind(cbind(yins_id[which(yins==yin)],ceng),weight), aldat)
  }
  aldat <- as.data.frame(list(case=aldat[,1],class=aldat[,2],weight=aldat[,3]))
  pics[[cnt]] <- ggplot(aldat, aes(case, weight, fill=class))+
    geom_bar(stat='identity',position=position_dodge(0.7),width=0.6)+
    theme_bw()+
    labs(title=timu[info])
  cnt <- cnt + 1
}

firstid <- c('B','D','E','F','G')
middleid <- c('H','I','J','K','L','W','X','AI','AO','AP','AR')
lastid <- c('AT','AU','AV')

first <- sapply(firstid, fromLetter)
middle <- sapply(middleid, fromLetter)
last <- sapply(lastid, fromLetter)
firstv <- PCA(dat[,first], draw = F)
middlev <- PCA(dat[,middle], draw = F)
lastv <- PCA(dat[,last], draw = F)

PCA(ndat)
